!--------------------------------
! Programa dedicado a calcular
! a integral de uma funcao
! de algumas formas diferentes.

! MODELO DE ARQUIVO DE ENTRADA AO FINAL
! DO CODIGO

! As seguintes regras para calcular a
! integral estão implementadas:
!
! Regra do trapezio:
! \int_-h^+h = (h/2)*(f_1+2*f_0+f_-1)
!
! Regra de Simpson:
! \int_-h^+h = (h/3)*(f_1+4*f_0+f_-1)
!
! Regra de Bode:
! \int_{x_0}^{x_4} = (2*h/45)*(7*f_0+32*f_1+12*f_2+32*f_3+7*f_4)


! Programa compilado e testado com:
! $ gfortran -o exerB exerB.f90
! $ ./exerA

! Forma de utilizacao auto-explicativa
! atravez de comentarios.

! _o_o_ oOo _o_o_ 23/10/2010


PROGRAM exerA
  IMPLICIT NONE
  INTEGER :: NN, i, apont, iintegB, iintegT, iintegS
  INTEGER, DIMENSION(20) :: N
  REAL (KIND=8), DIMENSION(20) :: H, pts
  REAL (KIND=8), DIMENSION(0:4) :: F
  REAL (KIND=8) :: integral, integB, integT, integS, bintegB=1E10, bintegT=1E10, bintegS=1E10

  ! Valor da integral calculado na mao
  integral =16./17.*exp(1./4.)*(1./4.*sin(1.)-cos(1.)) - 16./17.*(-cos(0.))
  write(*,*) "integral calculada na mao: ", integral

  ! Lendo os dados
  ! 1) Numero de intervalos
  OPEN(11,file='tab2_in.dat')
  read(11,*) NN

  ! 2) Intervalos em si
  i=1
  DO WHILE (i<=NN)
    read(11,*) N(i)
    H(i) = 1./N(i)
    i=i+1
  END DO

  open(unit=10,file="tab2_out.dat")
  write(10,*) "N(i)  |  H(i)  |  Regra do Trapezio  |  Regra de Simpson  |  Regra de Bode"
  ! Calculando as integrais aproximadas
  ! e os erros.
  i=1
  DO WHILE (i<=NN)
    apont=0
    integB=0
    integT=0
    integS=0

    ! Calculando a integral com um mesmo numero
    ! de divisoes.
    DO WHILE (apont<N(i))
      F(0)=exp( H(i)*(apont+0)/4. )*sin( H(i)*(apont+0))
      F(1)=exp( H(i)*(apont+1)/4. )*sin( H(i)*(apont+1))
      F(2)=exp( H(i)*(apont+2)/4. )*sin( H(i)*(apont+2))
      F(3)=exp( H(i)*(apont+3)/4. )*sin( H(i)*(apont+3))
      F(4)=exp( H(i)*(apont+4)/4. )*sin( H(i)*(apont+4))
      apont=apont+4

      ! Regra do Trapezio
      integT=integT+(F(0)+2*F(1)+2*F(2)+2*F(3)+F(4))

      ! Regra de Bode
      integB=integB+(7*F(0)+32*F(1)+12*F(2)+32*F(3)+7*F(4))

      ! Regra de Simpson
      integS=integS+(F(0)+4*F(1)+2*F(2)+4*F(3)+F(4))
    END DO

    ! Comparacoes para:
    ! 1) Registrar o erro de cada aproximacao
    ! 2) Tomar nota da melhor solucao
    integB = integB*(2*H(i)/45)  -integral
    IF (bintegB > integB) THEN
      bintegB=integB
      iintegB=i
    END IF

    integT=integT*(H(i)/2) -integral
    IF (bintegT > integT) THEN
      bintegT=integT
      iintegT=i
    END IF

    integS=integS*(H(i)/3) -integral
    IF (bintegS > integS) THEN
      bintegS=integS
      iintegS=i
    END IF

    write(10,*) N(i), H(i), integT, integS, integB
    i=i+1
  END DO

  write(*,*) "---------------------------------"
  write(*,*) "---------------------------------"
  write(*,*) "Para a regra do trapezio,"
  write(*,*) "o melhor valor de h eh", H(iintegT)
  write(*,*) "Pois o erro eh: ", bintegT

  write(*,*) "Para a regra de Simpson,"
  write(*,*) "o melhor valor de h eh", H(iintegS)
  write(*,*) "Pois o erro eh: ", bintegS

  write(*,*) "Para a regra de Bode,"
  write(*,*) "o melhor valor de h eh", H(iintegB)
  write(*,*) "Pois o erro eh: ", bintegB
  
end PROGRAM
  

!!!!!!!!!!!!!!!!!!
! MODELO DE ARQUIVO DE ENTRADA.
! DESCOMENTE e COPIE em tab2_in.dat
! 4
! 4
! 8
! 16
! 32
! 64
! 128
! 256
! 512
! 1024
! 2048
! 4096